#ifndef PCPS_LOCAL_ENERGY_HEADER
#define PCPS_LOCAL_ENERGY_HEADER

#include <src/pcps.h>
#include <src/input.h>
#include <src/operator.h>
#include <src/hamiltonian.h>
#include <src/wavefunction.h>
#include <src/coeff_value.h>

namespace pcps {

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   function object to compute the wavefunction's local energy
  ///
  ///    This class is templated on the scalar type (SCALAR), the particle type (PART),
  ///    and the reference type (REF).
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////

  template <class SCALAR, class PART, class REF> class LocalEnergyFunc {

    private:

      const pcps::Input & _userinp; ///< reference to the user input options
      const pcps::Hamiltonian<SCALAR> & _ham; ///< reference to the Hamiltonian
      const int _configs_per_site; ///< number of configurations on a single lattice site
      const int _nsites; ///< number of lattice sites
      const bool _qc; ///< whether we are using a quantum chemistry Hamiltonian
      const bool _use_enuc; ///< whether to include the nuclear (and core) energy
      const double _tei_thresh; ///< threshold for screening the two electron integrals
      std::vector<pcps::OpRL<SCALAR>      > _opsRL; ///< the Hamiltonian's raising/lowering operators
      std::vector<pcps::OpRLRL<SCALAR>    > _opsRLRL; ///< the Hamiltonian's double raising/lowering operators
      std::vector<pcps::OpNN<SCALAR>      > _opsNN; ///< the Hamiltonian's double number operators
      std::vector<pcps::OpN<SCALAR>       > _opsN; ///< the Hamiltonian's number operators
      std::vector<pcps::OpIdentity<SCALAR> > _opsIdentity; ///< the Hamiltonian's identity operators
      const pcps::action_value_object<SCALAR, pcps::OpRL<SCALAR>,       PART> _avoRL; ///< object for computing raising/lowering operator action
      const pcps::action_value_object<SCALAR, pcps::OpRLRL<SCALAR>,     PART> _avoRLRL; ///< object for computing double raising/lowering operator action
      const pcps::action_value_object<SCALAR, pcps::OpNN<SCALAR>,       PART> _avoNN; ///< object for computing double number operator action
      const pcps::action_value_object<SCALAR, pcps::OpN<SCALAR>,        PART> _avoN; ///< object for computing number operator action
      const pcps::action_value_object<SCALAR, pcps::OpIdentity<SCALAR>, PART> _avoIdentity; ///< object for computing identity operator action

      LocalEnergyFunc(); ///< default constructor is disabled
      LocalEnergyFunc(const LocalEnergyFunc &); ///< copy constructor is disabled
      LocalEnergyFunc & operator=(const LocalEnergyFunc &); ///< assignment operator is disabled

    public:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief   constructs the local energy function object
      ///
      /// \param   userinp   user input options
      /// \param   ham       the Hamiltonian
      /// \param   wfn       the wavefunction
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      LocalEnergyFunc(const pcps::Input & userinp,
                      const pcps::Hamiltonian<SCALAR> & ham,
                      const pcps::Wavefunction<SCALAR> & wfn)
        : _avoRL      (userinp.configs_per_site()),
          _avoRLRL    (userinp.configs_per_site()),
          _avoNN      (userinp.configs_per_site()),
          _avoN       (userinp.configs_per_site()),
          _avoIdentity(userinp.configs_per_site()),
          _userinp(userinp),
          _ham(ham),
          _configs_per_site(userinp.configs_per_site()),
          _nsites(userinp.nsites()),
          _qc( ham.is_qc() ),
          _use_enuc(!userinp.active_energy_only()),
          _tei_thresh(userinp.tei_screening())
      {
        // separate the hamiltonian's operators into their different types
        for (typename std::vector<typename pcps::OpBase<SCALAR>::const_ptr>::const_iterator op = ham.ops().begin(); op != ham.ops().end(); op++) {
          if ( typeid(**op) == typeid(pcps::OpRL<SCALAR>) ) {
            boost::shared_ptr<const pcps::OpRL<SCALAR> > derived_op_ptr
              = boost::dynamic_pointer_cast<const pcps::OpRL<SCALAR>, const pcps::OpBase<SCALAR> >(*op);
            _opsRL.push_back(*derived_op_ptr);
          } else if ( typeid(**op) == typeid(pcps::OpRLRL<SCALAR>) ) {
            boost::shared_ptr<const pcps::OpRLRL<SCALAR> > derived_op_ptr
              = boost::dynamic_pointer_cast<const pcps::OpRLRL<SCALAR>, const pcps::OpBase<SCALAR> >(*op);
            _opsRLRL.push_back(*derived_op_ptr);
          } else if ( typeid(**op) == typeid(pcps::OpNN<SCALAR>) ) {
            boost::shared_ptr<const pcps::OpNN<SCALAR> > derived_op_ptr
              = boost::dynamic_pointer_cast<const pcps::OpNN<SCALAR>, const pcps::OpBase<SCALAR> >(*op);
            _opsNN.push_back(*derived_op_ptr);
          } else if ( typeid(**op) == typeid(pcps::OpN<SCALAR>) ) {
            boost::shared_ptr<const pcps::OpN<SCALAR> > derived_op_ptr
              = boost::dynamic_pointer_cast<const pcps::OpN<SCALAR>, const pcps::OpBase<SCALAR> >(*op);
            _opsN.push_back(*derived_op_ptr);
          } else if ( typeid(**op) == typeid(pcps::OpIdentity<SCALAR>) ) {
            boost::shared_ptr<const pcps::OpIdentity<SCALAR> > derived_op_ptr
              = boost::dynamic_pointer_cast<const pcps::OpIdentity<SCALAR>, const pcps::OpBase<SCALAR> >(*op);
            _opsIdentity.push_back(*derived_op_ptr);
          } else
            throw pcps::Exception( (boost::format("unexpected operator type in LocalEnergyFunc constructor: \"%s\"") % typeid(**op).name()).str() );
        }
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief   computes the local energy for a given configuration
      ///
      /// \param[in]      config     the current lattice configuration, size _nsites
      /// \param[out]     workspace  space in which to manipulate the lattice configuration,
      ///                            size 2*_nsites
      /// \param[in,out]  cvf        object to compute wavefunction coefficient ratios, which must be
      ///                            initialized for the current configuration
      ///
      /// \return the local energy of the current configuration
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      SCALAR get_local_energy(const int * const config, int * const workspace, pcps::CoeffValueFunc<SCALAR, PART, REF> & cvf) const {

        // if we are dealing with a quantum chemistry Hamiltonian, return the corresponding local energy
        if (_qc) return this->get_qc_local_energy(config, workspace, cvf);

        // initialize the local energy
        SCALAR retval = pcps::zero<SCALAR>();

        // place the configuration in the workspace
        for (int i = 0; i < _nsites; i++)
          workspace[i] = config[i];

        // create an array to hold operators' sites' occupations
        int op_occs[4];

        // add each RL operator's contribution to the energy
        for (int i = 0; i < _opsRL.size(); i++) {

          // get a reference to the operator
          const pcps::OpRL<SCALAR> & op = _opsRL[i];

          // get the configuration of the ket and check whether it is valid
          bool valid = true;
          for (int k = 0; k < op.sites().size(); k++) {
            int & site_config = workspace[op.sites()[k]];
            site_config -= op.offsets()[k];
            op_occs[k] = site_config;
            if ( site_config < 0 || site_config >= _configs_per_site )
              valid = false;
          }

          // if the ket is a valid configuration, get the contribution to the energy
          // NOTE:  the sign from fermion operators is taken care of by the get_coeff_ratio function
          if (valid)
            retval += op.coeff() * _avoRL.get_value(op_occs, _configs_per_site)
                                 * cvf.get_coeff_ratio(2, &op.sites()[0], workspace, config);

          // revert back to the bra configuration
          for (int k = 0; k < op.sites().size(); k++)
            workspace[op.sites()[k]] += op.offsets()[k];

        }

        // add each RLRL operator's contribution to the energy
        for (int i = 0; i < _opsRLRL.size(); i++) {

          // get a reference to the operator
          const pcps::OpRLRL<SCALAR> & op = _opsRLRL[i];

          // get the configuration of the ket and check whether it is valid
          bool valid = true;
          for (int k = 0; k < op.sites().size(); k++) {
            int & site_config = workspace[op.sites()[k]];
            site_config -= op.offsets()[k];
            op_occs[k] = site_config;
            if ( site_config < 0 || site_config >= _configs_per_site )
              valid = false;
          }

          // if the ket is a valid configuration, get the contribution to the energy
          // NOTE:  the sign from fermion operators is taken care of by the get_coeff_ratio function
          if (valid)
            retval += op.coeff() * _avoRLRL.get_value(op_occs, _configs_per_site)
                                 * cvf.get_coeff_ratio(4, &op.sites()[0], workspace, config);

          // revert back to the bra configuration
          for (int k = 0; k < op.sites().size(); k++)
            workspace[op.sites()[k]] += op.offsets()[k];

        }

        // add each NN operator's contribution to the energy
        for (int i = 0; i < _opsNN.size(); i++) {

          // get a reference to the operator
          const pcps::OpNN<SCALAR> & op = _opsNN[i];

          // get the operator's sites' occupations
          op_occs[0] = config[op.sites()[0]];
          op_occs[1] = config[op.sites()[1]];

          // compute the energy contribution
          retval += op.coeff() * _avoNN.get_value(op_occs, _configs_per_site);

        }

        // add each N operator's contribution to the energy
        for (int i = 0; i < _opsN.size(); i++) {

          // get a reference to the operator
          const pcps::OpN<SCALAR> & op = _opsN[i];

          // get the operator's sites' occupations
          op_occs[0] = config[op.sites()[0]];

          // compute the energy contribution
          retval += op.coeff() * _avoN.get_value(op_occs, _configs_per_site);

        }

        // add each Identity operator's contribution to the energy
        for (int i = 0; i < _opsIdentity.size(); i++) {

          // get a reference to the operator
          const pcps::OpIdentity<SCALAR> & op = _opsIdentity[i];

          // compute the energy contribution
          retval += op.coeff() * _avoIdentity.get_value(op_occs, _configs_per_site);

        }

        // return the local energy
        return retval;

      }

    private:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief   computes the local energy for a given configuration for a quantum chemical Hamiltonian
      ///
      ///   There are 6 types of operators to consider.
      ///
      ///   From the one electron part of the Hamiltonian:
      ///
      ///     \f[  a^\dag_p a_p               \f]
      ///     \f[  a^\dag_p a_q               \f]
      ///
      ///   From the two electron part of the Hamiltonian:
      ///
      ///     \f[  a^\dag_p a_p a^\dag_q a_q  \f]
      ///     \f[  a^\dag_p a_p a^\dag_q a_r  \f]
      ///     \f[  a^\dag_p a_q a^\dag_r a_r  \f]
      ///     \f[  a^\dag_p a_q a^\dag_r a_s  \f]
      ///
      ///   For the one electron operators we have the alpha/beta combinations
      ///   aa, bb.
      ///
      ///   For the two electron operators we have the alpha/beta combinations
      ///   aaaa, aabb, bbaa, bbbb.
      ///
      ///   This all totals to 4 one electron cases and 16 two electron cases.
      ///
      /// \param[in]      config     the current lattice configuration, size _nsites
      /// \param[out]     workspace  space in which to manipulate the lattice configuration,
      ///                            size 2*_nsites
      /// \param[in,out]  cvf        object to compute wavefunction coefficient ratios, which must be
      ///                            initialized for the current configuration
      ///
      /// \return the local energy of the current configuration
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      SCALAR get_qc_local_energy(const int * const config, int * const workspace, pcps::CoeffValueFunc<SCALAR, PART, REF> & cvf) const {

        // get pointers to where we will store the singly and doubly excited configurations
        int * const sconf = workspace;
        int * const dconf = workspace + _nsites;

        // initialize the return value
        SCALAR retval = pcps::zero<SCALAR>();

        // compute aa, aaaa, aabb contributions
        retval += this->aaxx_contribution(0, config, sconf, dconf, cvf);

        // compute bb, bbbb, bbaa contributions
        retval += this->aaxx_contribution(1, config, sconf, dconf, cvf);

        // add the nuclear repulsion energy
        if (_use_enuc)
          retval += pcps::unity<SCALAR>() * _ham.enuc();

        // return the local energy
        return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief   computes the local energy contributions when first excitation is of the given spin
      ///
      /// \param[in]      offset     sets the spin of the first excitation operator.  0 = alpha, 1 = beta
      /// \param[in]      config     the current lattice configuration, size _nsites
      /// \param[out]     sconf      configuration after the 1st excitation
      /// \param[out]     dconf      configuration after the 2nd excitation
      /// \param[in,out]  cvf        object to compute wavefunction coefficient ratios, which must be
      ///                            initialized for the current configuration
      ///
      /// \return the (aa, aaaa, and aabb) or (bb, bbbb, and bbaa) contributions to the local energy
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      SCALAR aaxx_contribution(const int offset, const int * const config, int * const sconf, int * const dconf,
                               pcps::CoeffValueFunc<SCALAR, PART, REF> & cvf) const {

        // get the number of quantum chemistry orbitals
        const int n = _nsites / 2;

        // create an array for the sites being affected by an excitation
        int affected[2];

        // initialize the singly excited configuration
        for (int i = 0; i < _nsites; i++)
          sconf[i] = config[i];

        // initialize the return value
        SCALAR retval = pcps::zero<SCALAR>();

        // compute one-electron contribution
        for (int p = 0; p < n; p++)
        for (int q = 0; q < n; q++) {

          // number operators are treated specially
          if ( p == q ) {
            retval += ( sconf[2*p+offset] == 1 ? _ham.get_mei(p,q) : pcps::zero<SCALAR>() );
            continue;
          }

          // create the excitation, skipping excitations that destroy the configuration
          if ( sconf[2*p+offset] == 1 && sconf[2*q+offset] == 0 ) {
            sconf[2*p+offset] = 0;
            sconf[2*q+offset] = 1;
          } else {
            continue;
          }

          // compute the excitation's contribution
          affected[0] = 2*p+offset;
          affected[1] = 2*q+offset;
          retval += _ham.get_mei(p,q) * cvf.get_coeff_ratio(2, affected, sconf, config);

          // remove the excitation
          sconf[2*p+offset] = 1;
          sconf[2*q+offset] = 0;

        }

        // ensure the length of the vector of non-zero pairs is sane
        if (_ham.nz_pairs().size() % 2 != 0)
          throw pcps::Exception("number of non-zero pair indices must be even in pcps::LocalEnergyFunc::aaxx_contribution");

        long int tei_count = 0;

        // compute two-electron contribution by looping over non-zero excitation pairs
        for (std::vector<int>::const_iterator outer_it = _ham.nz_pairs().begin(); outer_it != _ham.nz_pairs().end(); ) {

          // get the orbital indices for the first excitation
          const int p = *outer_it++;
          const int q = *outer_it++;

          // create the excitation, skipping excitations that destroy the configuration
          if (p == q) {
            if ( sconf[2*p+offset] == 0 )
              continue;
          } else if ( sconf[2*p+offset] == 1 && sconf[2*q+offset] == 0 ) {
            sconf[2*p+offset] = 0;
            sconf[2*q+offset] = 1;
          } else {
            continue;
          }

          // get the wavefunction ratio between the base and first excited configuration
          affected[0] = 2*p+offset;
          affected[1] = 2*q+offset;
          const SCALAR outer_ratio = ( p == q ? pcps::unity<SCALAR>() : cvf.get_coeff_ratio(2, affected, sconf, config) );

          // initialize the doubly excited configuration
          for (int i = 0; i < _nsites; i++)
            dconf[i] = sconf[i];

          // if the wavefunction is non-zero at the configuration of the first excitation,
          // update the cvf object and compute ratios for the second excitation
          if ( std::abs(outer_ratio) > 1.0e-9 ) {

            // update the coefficient value function object for the excited configuration
            affected[0] = 2*p+offset;
            affected[1] = 2*q+offset;
            if (p != q)
              cvf.update(2, affected, sconf, config);

            // loop over the second excitations
            for (std::vector<int>::const_iterator inner_it = _ham.nz_pairs().begin(); inner_it != _ham.nz_pairs().end(); ) {

              // get the orbital indices for the second excitation
              const int r = *inner_it++;
              const int s = *inner_it++;

              // get the two-electron integral
              const SCALAR tei = _ham.get_tei(p,q,r,s);

              // skip integrals that are below the screening threshold
              if ( std::abs(tei) < _tei_thresh ) continue;

              // get the contribution when the second excitation is alpha
              retval += 0.5 * outer_ratio * tei * this->inner_ratio(0, r, s, affected, sconf, dconf, cvf, tei_count);

              // get the contribution when the second excitation is beta
              retval += 0.5 * outer_ratio * tei * this->inner_ratio(1, r, s, affected, sconf, dconf, cvf, tei_count);

            }

            // update the coefficient value function object for the de-excited configuration
            affected[0] = 2*p+offset;
            affected[1] = 2*q+offset;
            if (p != q)
              cvf.update(2, affected, config, sconf);

          // if the wavefunction is zero at the configuration for the first excitation,
          // directly compute the wavefunction ratios for each of the second excitation configurations
          } else {

            // get the wavefunction value for the initial configuration in mantissa/exponent format
            pcps::MantExp<SCALAR> initial_value = cvf.get_wfn_value(config);

            // include the sign caused by the first excitation operator
            initial_value *= this->qc_op_sign(offset, p, q, config);

            SCALAR contribution = pcps::zero<SCALAR>();

            // loop over the second excitations
            for (std::vector<int>::const_iterator inner_it = _ham.nz_pairs().begin(); inner_it != _ham.nz_pairs().end(); ) {

              // get the orbital indices for the second excitation
              const int r = *inner_it++;
              const int s = *inner_it++;

              // get the two-electron integral
              const SCALAR tei = _ham.get_tei(p,q,r,s);

              // skip integrals that are below the screening threshold
              if ( std::abs(tei) < _tei_thresh ) continue;

              // number operators
              if (r == s) {

                // get the contribution when the second excitation is alpha
                contribution += 0.5 * tei * outer_ratio * double(dconf[2*r+0]);

                // get the contribution when the second excitation is beta
                contribution += 0.5 * tei * outer_ratio * double(dconf[2*r+1]);

              // excitation operators
              } else {

                // get the contribution when the second excitation is alpha
                contribution += 0.5 * tei * this->direct_ratio(0, r, s, initial_value, dconf, cvf, tei_count);

                // get the contribution when the second excitation is beta
                contribution += 0.5 * tei * this->direct_ratio(1, r, s, initial_value, dconf, cvf, tei_count);

              }

            }

            retval += contribution;

            //pthread_mutex_lock(pcps::alloc_mutex());
            //std::ofstream f("special_doubles.txt", std::ios_base::app);
            //for (int i = 0; i < _nsites; i++)
            //  f << config[i];
            //f << boost::format(" %2i %2i %s") % p % q % pcps::format_number("%12.4e", contribution) << std::endl;
            //f.close();
            //pthread_mutex_unlock(pcps::alloc_mutex());

          }

          // remove the excitation
          if (p != q) {
            sconf[2*p+offset] = 1;
            sconf[2*q+offset] = 0;
          }

        }

        //std::cout << "handled " << tei_count << " four index alpha operators in aaxx_contribution" << std::endl;

        // return the result
        return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief   computes the ratio between the second and first excitations' wavefunction values
      ///
      /// \param[in]      offset     sets the spin of the second excitation operator.  0 = alpha, 1 = beta
      /// \param[in]      r          1st orbital index of the second excitation
      /// \param[in]      s          2nd orbital index of the second excitation
      /// \param[out]     affected   array in which to store the sites touched by the second excitation
      /// \param[in]      sconf      configuration after the 1st excitation
      /// \param[in,out]  dconf      on entry and on exit this array holds the configuration after the
      ///                            first excitation, but during the function excecution it is changed
      ///                            to the configuration after the second excitation
      /// \param[in,out]  cvf        object to compute wavefunction coefficient ratios,
      ///                            which must be initialized for the configuration
      ///                            resulting from the first exciataion
      ///
      /// \return the ratio between the second and first excitations' wavefunction values
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      SCALAR inner_ratio(const int offset, const int r, const int s,
                         int * const affected, const int * const sconf, int * const dconf,
                         pcps::CoeffValueFunc<SCALAR, PART, REF> & cvf,
                         long int & tei_count) const {

        // number operators are treated specially
        if ( r == s )
          return ( dconf[2*r+offset] == 1 ? pcps::unity<SCALAR>() : pcps::zero<SCALAR>() );

        // create the excitation, skipping excitations that destroy the configuration
        if ( dconf[2*r+offset] == 1 && dconf[2*s+offset] == 0 ) {
          dconf[2*r+offset] = 0;
          dconf[2*s+offset] = 1;
        } else {
          return pcps::zero<SCALAR>();
        }

        tei_count += (long int)offset;

        // get the wavefunction ratio between the second and fisrt excited configurations
        affected[0] = 2*r+offset;
        affected[1] = 2*s+offset;
        const SCALAR retval = cvf.get_coeff_ratio(2, affected, dconf, sconf);

        // remove the excitation
        dconf[2*r+offset] = 1;
        dconf[2*s+offset] = 0;

        // return the ratio
        return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief   computes the wavefunction ratio between the doubly excited and initial configurations
      ///
      /// \param[in]      offset     sets the spin of the second excitation operator.  0 = alpha, 1 = beta
      /// \param[in]      r          1st orbital index of the second excitation
      /// \param[in]      s          2nd orbital index of the second excitation
      /// \param[in]      iv         Wavefunction value for the initial configuration
      /// \param[in,out]  dconf      on entry and on exit this array holds the configuration after the
      ///                            first excitation, but during the function excecution it is changed
      ///                            to the configuration after the second excitation
      /// \param[in,out]  cvf        object to compute wavefunction coefficient values
      ///
      /// \return the wavefunction ratio between the doubly excited and initial configurations
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      SCALAR direct_ratio(const int offset,
                          const int r,
                          const int s,
                          const pcps::MantExp<SCALAR> & iv,
                          int * const dconf,
                          pcps::CoeffValueFunc<SCALAR, PART, REF> & cvf,
                          long int & tei_count) const
      {

        // number operators give zero, as this function is only called when the wavefunction coefficient
        // for the first excited configuration is negligible
        if ( r == s )
          return pcps::zero<SCALAR>();

        // create the excitation, skipping excitations that destroy the configuration
        if ( dconf[2*r+offset] == 1 && dconf[2*s+offset] == 0 ) {
          dconf[2*r+offset] = 0;
          dconf[2*s+offset] = 1;
        } else {
          return pcps::zero<SCALAR>();
        }

        tei_count += (long int)offset;

        // get the ratio of wavefunction values for the doubly excited and initial configurations in mantissa/exponent
        pcps::MantExp<SCALAR> ratio = cvf.get_wfn_value(dconf) / iv;

        // remove the excitation
        dconf[2*r+offset] = 1;
        dconf[2*s+offset] = 0;

        // include the sign caused by the second excitation operator
        ratio *= this->qc_op_sign(offset, r, s, dconf);

        // return the ratio as a regular number
        return ratio.get();

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief   computes the sign resulting from the operation of an excitation operator
      ///
      /// \param[in]      offset     sets the spin of the excitation operator.  0 = alpha, 1 = beta
      /// \param[in]      r          1st orbital index of the excitation
      /// \param[in]      s          2nd orbital index of the excitation
      /// \param[in]      config     the configuration of the bra state
      ///
      /// \return the sign produced by operating the excitation operator on the bra
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      double qc_op_sign(const int offset, const int r, const int s, const int * const config) const {

        // get max/min markers denoting the range of sites between r and s
        const int i_min = 2 * std::min(r,s) + offset;
        const int i_max = 2 * std::max(r,s) + offset;

        // initialize the sign
        double retval = 1.0; //( r <= s ? 1.0 : -1.0 );

        // For each occupied site between the excitation operator's sites, add a sign of -1.
        // The behavior is different for different references, as they rely on different definitions
        // of the normal ordering of alpha and beta sites.
        if ( typeid(REF) == typeid(pcps::Ref_Uniform) || typeid(REF) == typeid(pcps::Ref_Pfaffian) ) {

          // Uniform and Pfaffian references work with an ababab ordering
          for (int i = i_min+1; i < i_max; i++)
            retval *= 1 - 2 * config[i];

        } else {

          // RHF, AGP, and Pairing references work with an aaabbb ordering
          for (int i = i_min+2; i < i_max; i+=2)
            retval *= 1 - 2 * config[i];

        }

        // return the total sign
        return retval;

      }

  };

} // end namespace pcps

#endif
